

         g_up_xx = auxfields.u_aux[G_UP_XX]GRD_IND;
         g_up_yy = auxfields.u_aux[G_UP_YY]GRD_IND;
         g_up_xy = auxfields.u_aux[G_UP_XY]GRD_IND;

	 detg_1_2 = 1./sqrt(g_up_xx*g_up_yy - g_up_xy*g_up_xy);

#ifdef PROJECTION



       m_x = (*fields_ptr).u[M_A]GRD_IND;
       m_y = (*fields_ptr).u[M_B]GRD_IND;
       m_z = (*fields_ptr).u[M_C]GRD_IND;



       m_x_t =  (dfields.du[DD_X_M_A]GRD_IND +dfields.du[DD_Y_M_A]GRD_IND)/detg_1_2 ;
       m_y_t =  (dfields.du[DD_X_M_B]GRD_IND +dfields.du[DD_Y_M_B]GRD_IND)/detg_1_2 ;
       m_z_t =  (dfields.du[DD_X_M_C]GRD_IND +dfields.du[DD_Y_M_C]GRD_IND)/detg_1_2 ;

       m_t_dot_m = m_x_t*m_x + m_y_t*m_y + m_z_t*m_z;
       norm_m = m_x*m_x + m_y*m_y + m_z*m_z;

       m_x_t = m_x_t - m_x*m_t_dot_m/norm_m;
       m_y_t = m_y_t - m_y*m_t_dot_m/norm_m;
       m_z_t = m_z_t - m_z*m_t_dot_m/norm_m;




#else


       d2_x =    g_up_xx*dfields.du[D_X_M_A]GRD_IND*dfields.du[D_X_M_A]GRD_IND +
              2.*g_up_xy*dfields.du[D_X_M_A]GRD_IND*dfields.du[D_Y_M_A]GRD_IND +
                 g_up_yy*dfields.du[D_Y_M_A]GRD_IND*dfields.du[D_Y_M_A]GRD_IND;
       d2_y =    g_up_xx*dfields.du[D_X_M_B]GRD_IND*dfields.du[D_X_M_B]GRD_IND +
              2.*g_up_xy*dfields.du[D_X_M_B]GRD_IND*dfields.du[D_Y_M_B]GRD_IND +
                 g_up_yy*dfields.du[D_Y_M_B]GRD_IND*dfields.du[D_Y_M_B]GRD_IND;
       d2_z =    g_up_xx*dfields.du[D_X_M_C]GRD_IND*dfields.du[D_X_M_C]GRD_IND +
              2.*g_up_xy*dfields.du[D_X_M_C]GRD_IND*dfields.du[D_Y_M_C]GRD_IND +
                 g_up_yy*dfields.du[D_Y_M_C]GRD_IND*dfields.du[D_Y_M_C]GRD_IND;

       m_x = (*fields_ptr).u[M_A]GRD_IND;
       m_y = (*fields_ptr).u[M_B]GRD_IND;
       m_z = (*fields_ptr).u[M_C]GRD_IND;

       md_xm = m_x*dfields.du[D_X_M_A]GRD_IND + m_y*dfields.du[D_X_M_B]GRD_IND + m_z*dfields.du[D_X_M_C]GRD_IND; 

       md_ym = m_x*dfields.du[D_Y_M_A]GRD_IND + m_y*dfields.du[D_Y_M_B]GRD_IND + m_z*dfields.du[D_Y_M_C]GRD_IND; 

       mmd2 = g_up_xx*md_xm*md_xm + 2.*g_up_xy*md_xm*md_ym + g_up_yy*md_ym*md_ym;

       d2 = -(d2_x + d2_y + d2_z - mmd2) - mm*(1. - m_x*m_x - m_y*m_y - m_z*m_z);

       m_x_t =  (dfields.du[DD_X_M_A]GRD_IND +dfields.du[DD_Y_M_A]GRD_IND)/detg_1_2 - m_x*d2;
       m_y_t =  (dfields.du[DD_X_M_B]GRD_IND +dfields.du[DD_Y_M_B]GRD_IND)/detg_1_2 - m_y*d2;
       m_z_t =  (dfields.du[DD_X_M_C]GRD_IND +dfields.du[DD_Y_M_C]GRD_IND)/detg_1_2 - m_z*d2;

#endif


#ifdef DISSIPATION

  m_x_t = m_x_t + sigma*dfields.du[DISS_M_A]GRD_IND;
  m_y_t = m_y_t + sigma*dfields.du[DISS_M_B]GRD_IND;
  m_z_t = m_z_t + sigma*dfields.du[DISS_M_C]GRD_IND;

#endif

/* printf("I am inside eq_eq"); */






// just wrong garbage


//#ifdef PROJECTION


	      //x =  -1.0 + (FLOAT)(grid_ind1)*one_dN1;
	      //y =  -1.0 + (FLOAT)(grid_ind2)*one_dN2;
                    
                    //Normal_cub(x,y,nn_x,nn_y,nn_z,(*grid_ptr).grid);
 


       //m_x_t =  (dfields.du[DD_X_M_A]GRD_IND +dfields.du[DD_Y_M_A]GRD_IND)/detg_1_2 ;
       //m_y_t =  (dfields.du[DD_X_M_B]GRD_IND +dfields.du[DD_Y_M_B]GRD_IND)/detg_1_2 ;
       //m_z_t =  (dfields.du[DD_X_M_C]GRD_IND +dfields.du[DD_Y_M_C]GRD_IND)/detg_1_2 ;

       //m_t_nn = m_x_t*nn_x + m_y_t*nn_y + m_z_t*nn_z;

       //m_x_t = m_x_t - nn_x*m_t_nn;
       //m_y_t = m_y_t - nn_y*m_t_nn;
       //m_z_t = m_z_t - nn_z*m_t_nn;




//#else
